/*
* Copyright (c) 2008-2011, Thomas Handorf, Nils Christian
* 
* This file is part of libScopes.
* libScopes is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* libScopes is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU Lesser General Public License
* along with libScopes.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "net.h"

Net::Net(){
	nstr=new NetStruct();
	nstr->refCount++;
	reaMask=new VC();
	
}
Net::Net(NetStruct  &str2){
	nstr=&str2;
	nstr->refCount++;
	reaMask=new VC(nstr->numR);
}
Net::Net(size_type nr, size_type nc){
	nstr=new NetStruct(nr,nc);
	nstr->refCount++;
	reaMask=new VC(nr);
}
Net::Net(const Net& net2){
	nstr=net2.nstr;
	nstr->refCount++;
	reaMask=new VC((*net2.reaMask));
}
/*Net::Net(string &st){
	size_type pos=0;
	nstr=new NetStruct(st,pos);
	nstr->refCount++;
	reaMask=new VC(st,pos);
}
Net::Net(string &st,size_type pos){
	nstr=new NetStruct(st,pos);
	nstr->refCount++;
	reaMask=new VC(st,pos);
}*/
Net::~Net(){
	nstr->refCount--;
	if (nstr->refCount==0) delete nstr;
	delete reaMask;
   reaMask=NULL;
#ifdef DEBUG 
   printf("Net::destructor\n");fflush(stdout);
#endif
}
Net& Net::operator= (const Net& net2){
   nstr->refCount--;
   if (nstr->refCount==0) delete nstr;
   delete reaMask;
   nstr=net2.nstr;
   nstr->refCount++;
   reaMask=new VC((*net2.reaMask));
   return *this;
}
void Net::onChange(){ // deferred copy of Netstruct; this is only invoked when the structure is changed after more than one net references the same netstruct
   //note: changes to reaMask should not trigger onChange()!!!!
   if (nstr->refCount>1){ // there are other net objects on this nstr -> need to copy; 
      nstr->refCount--; // that's us
      nstr=new NetStruct(*nstr); // copy whole NetStruct;
   }
}

/*void Net::initFromString(string &st,size_type pos){
	if (nstr){
		nstr->refCount--;
		if (nstr->refCount==0) delete nstr;
	}
	if (reaMask){
		delete reaMask;
	}
	nstr=new NetStruct(st,pos);
	nstr->refCount++;
	reaMask=new VC(st,pos);
}*/
/*!
	\brief sets the network size.
	\param nr number of reactions
	\param nc number of compounds
	
	Resizes the underlying array accordingly, thus avoiding continious reallocation if the reactions are added one by one.
*/
void Net::resize(size_type nr, size_type nc){
	if ((nc!=nstr->numC) || (nr!=nstr->numR)){
      onChange();
		nstr->resize(nr,nc);
      reaMask->resize(nr);
	}
}

/*!
	\brief add a new reaction to the network.
	\param subs array of substrate indices
	\param prds array of product indices
	\param rev reversibility, 0=reversible, 1=irreversible, -1=irreversible but in prds->subs direction.
	\param subsst stoichiometries for the substrates (same order). Omit if all stoich. equal 1
	\param prdsst stoichiometries for the products (same order). Omit if all stoich. equal 1
	\param cats catalysts of reaction
	\param inhs inhibitors of the reaction
	\param inhd inhibited substrates of the reaction (competitive)
	\param prps depletion propagators of the reactions
	
	places a new reaction at the next position after the lastly added reaction. if a reaction with that index existed it is removed.
	Note: all inhd's must also be substrates and all prps must also be products. If prps is given, substrate inhibition only occurs if at least one of the propagators is depleted
	
	The content of the arrays is copied. They can be deleted afterwards.
*/
size_type Net::addRea(VS* subs, VS* prds,char rev,VF* subsst,VF* prdsst,VS* cats,VS* inhs,VS* inhd,VS* prps){
	return replaceRea(nstr->nextToAdd,subs,prds,rev,subsst,prdsst,cats,inhs,inhd,prps);
}

/*!
	\brief place a new reaction into the network
	\param reaidx position where to place
	\param subs array of substrate indices
	\param prds array of product indices
	\param rev reversibility, 0=reversible, 1=irreversible, -1=irreversible but in prds->subs direction.
	\param subsst stoichiometries for the substrates (same order). Omit if all stoich. equal 1
	\param prdsst stoichiometries for the products (same order). Omit if all stoich. equal 1
	\param cats catalysts of reaction
	\param inhs inhibitors of the reaction
	\param inhd inhibited substrates of the reaction (competitive)
	\param prps depletion propagators of the reactions
	
	places a new reaction at position reaidx. if a reaction with that index existed it is removed.
	Note: all inhd's must also be substrates and all prps must also be products. If prps is given, substrate inhibition only occurs if at least one of the propagators is depleted
	
	The content of the arrays is copied. They can be deleted afterwards.
*/
size_type Net::replaceRea(size_type reaidx,VS* subs, VS* prds, char rev,VF* subsst,VF* prdsst,VS* cats,VS* inhs,VS* inhd,VS* prps){
	size_type nsubs=subs->size();
	size_type nprds=prds->size();
	size_type i;
	bool remsst=false,rempst=false;
   onChange();
	if (!subsst){
		subsst=new VF();
		remsst=true; //remove generated array later
	}
	if (!prdsst){
		prdsst=new VF();
		rempst=true; //remove generated array later
	}
	if (nsubs!=subsst->size()){
		subsst->assign(nsubs,1); // create substrate stoichiometrie =1  if not provided as param
	}
	if (nprds!=prdsst->size()){
      prdsst->assign(nprds,1); // create product stoichiometrie =1 if not provided as param
	}
	size_type maxc=std::max(subs->max(),prds->max()); // search for maximum compund index
   if (inhs) maxc=std::max(maxc,inhs->max()); // search inhibitors for maximum compound index
   if (cats) maxc=std::max(maxc,cats->max()); // same for catalysts
   if (reaidx>=nstr->numR){ 
		if (maxc<nstr->numC) maxc=nstr->numC-1;
		resize(reaidx+1,maxc+1);
	} else {
		if (maxc>=nstr->numC) resize(nstr->numR,maxc+1);
	}
	//we del reaction first to update CSUC,CPRE
	delRea(reaidx);
	(*nstr->IN)[reaidx]->assign(*subs);
   (*nstr->OUT)[reaidx]->assign(*prds);
	if (inhs) (*nstr->INH)[reaidx]->assign(*inhs);
	if (inhd) (*nstr->CMPI)[reaidx]->assign(*inhd);
   if (cats) (*nstr->CAT)[reaidx]->assign(*cats);
   if (prps) (*nstr->PRP)[reaidx]->assign(*prps);
   (*nstr->INST)[reaidx]->assign(*subsst);
	(*nstr->OUTST)[reaidx]->assign(*prdsst);
	for (i=0;i<nsubs;i++){
		(*nstr->CSUC)[(*subs)[i]]->push_back(reaidx);
	}
	if (inhs){
		for (i=0;i<inhs->size();i++){
			(*nstr->CSUC)[(*inhs)[i]]->push_back(reaidx);
		}
	}
	if (cats){
      for (i=0;i<cats->size();i++){
         (*nstr->CSUC)[(*cats)[i]]->push_back(reaidx);
      }
   }
   for (i=0;i<nprds;i++){
		(*nstr->CPRE)[(*prds)[i]]->push_back(reaidx);
	}
	(*nstr->REV)[reaidx]=rev;
	(*reaMask)[reaidx]=1;
   if (reaidx>=nstr->nextToAdd) nstr->nextToAdd=reaidx+1;
	if (remsst) delete subsst;
	if (rempst) delete prdsst;
	return reaidx;
}
int Net::hasRea(VS* subs, VS* prds, char rev,VF* subsst,VF* prdsst){
	return 0; //not implemented yet
}

void Net::setSubstrates(size_type reaidx,VS* subs,VF* subsst){
   size_type nsubs=subs->size();
   size_type i;
   int pos;
   bool delst=false;
   onChange();
   if (!subsst){
      subsst=new VF();
      delst=true; //remove generated array later
   }
   if (nsubs!=subsst->size()){
      subsst->assign(nsubs,1); // create substrate stoichiometrie =1  if not provided as param
   }
   size_type maxc=subs->max();
   if (maxc>nstr->numC-1 || reaidx>nstr->numR-1){ // make sure network is large enoguh
      resize(std::max(reaidx+1,nstr->numR),std::max(maxc+1,nstr->numC));
   }
   for (i=0;i<(*nstr->IN)[reaidx]->size();i++){ // remove back links from old substrates
      pos=(*nstr->CSUC)[(*(*nstr->IN)[reaidx])[i]]->find(reaidx);
      if (pos>=0){
         (*nstr->CSUC)[(*(*nstr->IN)[reaidx])[i]]->erase(pos,1);
      }
   }
   for (i=0;i<nsubs;i++){
      (*nstr->CSUC)[(*subs)[i]]->push_back(reaidx);
   }
   (*nstr->IN)[reaidx]->assign(*subs);
   (*nstr->INST)[reaidx]->assign(*subsst);
   if (delst) delete subsst;
}
void Net::setProducts(size_type reaidx,VS* prds,VF* prdsst){
   size_type nprds=prds->size();
   size_type i;
   int pos;
   bool delst=false;
   onChange();
   if (!prdsst){
      prdsst=new VF();
      delst=true; //remove generated array later
   }
   if (nprds!=prdsst->size()){
      prdsst->assign(nprds,1); // create products stoichiometry =1  if not provided as param
   }
   size_type maxc=prds->max();
   if (maxc>nstr->numC-1 || reaidx>nstr->numR-1){ // make sure network is large enoguh
      resize(std::max(reaidx+1,nstr->numR),std::max(maxc+1,nstr->numC));
   }
   for (i=0;i<(*nstr->OUT)[reaidx]->size();i++){ // remove back links from old products
      pos=(*nstr->CPRE)[(*(*nstr->OUT)[reaidx])[i]]->find(reaidx);
      if (pos>=0){
         (*nstr->CPRE)[(*(*nstr->OUT)[reaidx])[i]]->erase(pos,1);
      }
   }
   for (i=0;i<nprds;i++){
      (*nstr->CPRE)[(*prds)[i]]->push_back(reaidx);
   }
   (*nstr->OUT)[reaidx]->assign(*prds);
   (*nstr->OUTST)[reaidx]->assign(*prdsst);
   if (delst) delete prdsst;
}
void Net::setInhibitors(size_type reaidx,VS* inh){
   size_type ninh=inh->size();
   size_type i;
   int pos;
   size_type maxc=inh->max();
   onChange();
   if (maxc>nstr->numC-1 || reaidx>nstr->numR-1){ // make sure network is large enoguh
      resize(std::max(reaidx+1,nstr->numR),std::max(maxc+1,nstr->numC));
   }
   for (i=0;i<(*nstr->INH)[reaidx]->size();i++){ // remove back links from old inhibitors
      pos=(*nstr->CSUC)[(*(*nstr->INH)[reaidx])[i]]->find(reaidx);
      if (pos>=0){
         (*nstr->CSUC)[(*(*nstr->INH)[reaidx])[i]]->erase(pos,1);
      }
   }
   for (i=0;i<ninh;i++){
      (*nstr->CSUC)[(*inh)[i]]->push_back(reaidx);
   }
   (*nstr->INH)[reaidx]->assign(*inh);
}
void Net::setCatalysts(size_type reaidx,VS* cat){
   size_type ncat=cat->size();
   size_type i;
   int pos;
   size_type maxc=cat->max();
   onChange();
   if (maxc>nstr->numC-1 || reaidx>nstr->numR-1){ // make sure network is large enoguh
      resize(std::max(reaidx+1,nstr->numR),std::max(maxc+1,nstr->numC));
   }
   for (i=0;i<(*nstr->CAT)[reaidx]->size();i++){ // remove back links from old catalysts
      pos=(*nstr->CSUC)[(*(*nstr->CAT)[reaidx])[i]]->find(reaidx);
      if (pos>=0){
         (*nstr->CSUC)[(*(*nstr->CAT)[reaidx])[i]]->erase(pos,1);
      }
   }
   for (i=0;i<ncat;i++){
      (*nstr->CSUC)[(*cat)[i]]->push_back(reaidx);
   }
   (*nstr->CAT)[reaidx]->assign(*cat);
   
}
void Net::setInhibited(size_type reaidx,VS* inhd){
   size_type maxc=inhd->max();
   onChange();
   if (maxc>nstr->numC-1 || reaidx>nstr->numR-1){ // make sure network is large enoguh
      resize(std::max(reaidx+1,nstr->numR),std::max(maxc+1,nstr->numC));
   }
   //no backlinks for CMPI
   (*nstr->CMPI)[reaidx]->assign(*inhd);
}  
void Net::setPropagators(size_type reaidx,VS* prp){
   size_type maxc=prp->max();
   onChange();
   if (maxc>nstr->numC-1 || reaidx>nstr->numR-1){ // make sure network is large enoguh
      resize(std::max(reaidx+1,nstr->numR),std::max(maxc+1,nstr->numC));
   }
   //no backlinks for PRP
   (*nstr->PRP)[reaidx]->assign(*prp);
}
void Net::setMembers(size_type comidx,VS* idxs){
   size_type max=idxs->max();
   onChange();
   if (comidx>max) max=comidx ;
   if (max>=nstr->numC) resize(nstr->numR,max+1);
   (*nstr->FAM)[comidx]->assign(*idxs);
   nstr->hasFam=1;
}
void Net::setChannels(size_type comidx,VS* idxs){
   onChange();
   if (comidx>=nstr->numC) resize(nstr->numR,comidx+1);
   (*nstr->FCN)[comidx]->assign(*idxs);
}
void Net::setName(size_type comidx,const std::string& name){
   onChange();
   if (comidx>=nstr->numC) resize(nstr->numR,comidx+1);
   (*nstr->names)[comidx]=name;
}
/*!
	\brief delete reaction at position idx
	\param idx index of reaction to delete
*/
void Net::delRea(size_type idx){
	size_type i;
	int pos;
   onChange();
	for (i=0;i<(*nstr->IN)[idx]->size();i++){
		pos=(*nstr->CSUC)[(*(*nstr->IN)[idx])[i]]->find(idx);
		if (pos>=0){
			(*nstr->CSUC)[(*(*nstr->IN)[idx])[i]]->erase(pos,1);
		}
	}
	for (i=0;i<(*nstr->INH)[idx]->size();i++){
		pos=(*nstr->CSUC)[(*(*nstr->INH)[idx])[i]]->find(idx);
		if (pos>=0){
			(*nstr->CSUC)[(*(*nstr->INH)[idx])[i]]->erase(pos,1);
		}
	}
	for (i=0;i<(*nstr->CAT)[idx]->size();i++){
      pos=(*nstr->CSUC)[(*(*nstr->CAT)[idx])[i]]->find(idx);
      if (pos>=0){
         (*nstr->CSUC)[(*(*nstr->CAT)[idx])[i]]->erase(pos,1);
      }
   }
   for (i=0;i<(*nstr->OUT)[idx]->size();i++){
		pos=(*nstr->CPRE)[(*(*nstr->OUT)[idx])[i]]->find(idx);
		if (pos>=0){
			(*nstr->CPRE)[(*(*nstr->OUT)[idx])[i]]->erase(pos,1);
		}
	}
	(*nstr->IN)[idx]->clear();
	(*nstr->OUT)[idx]->clear();
	(*nstr->INH)[idx]->clear();
	(*nstr->CMPI)[idx]->clear();
   (*nstr->CAT)[idx]->clear();
   (*nstr->INST)[idx]->clear();
	(*nstr->OUTST)[idx]->clear();
	(*reaMask)[idx]=0;
}
/*!
	\param reaidx which reaction
	\param comidx which compound
	
	returns the corresponding entry of the stoichiometric matrix
*/
float Net::at(size_type reaidx,size_type comidx){ // modify return value to double later
	int pos;
	float st=0; //modify to double later
	pos=(*nstr->IN)[reaidx]->find(comidx);
	if (pos>=0) {st=-(*(*nstr->INST)[reaidx])[pos];}
	pos=(*nstr->OUT)[reaidx]->find(comidx);
	if (pos>=0) {st+=(*(*nstr->OUTST)[reaidx])[pos];}
	return st;
}

/*!
	\param comidx which compound
	
	Note: this returns a newly generated array. Delete youself!
*/
VS* Net::getNeighbors(size_type comidx){
//	printf("comidx %li\n",comidx);fflush(stdout);
//	if (((*nstr->CSUC)[comidx]==NULL) ||((*nstr->CPRE)[comidx]==NULL)) printf("error null ptr\n");fflush(stdout);
	VS* ret=new VS(*(*nstr->CPRE)[comidx]);
   ret->put(ret->size(),*(*nstr->CSUC)[comidx]);
	return ret;
}
/*!
	\param reaidx which compound
	
	Note: this returns a newly generated array. Delete youself!
*/
VS* Net::getMetabolites(size_type reaidx){
//	printf("comidx %li\n",comidx);fflush(stdout);
//	if (((*nstr->CSUC)[comidx]==NULL) ||((*nstr->CPRE)[comidx]==NULL)) printf("error null ptr\n");fflush(stdout);
	VS* ret=new VS(*(*nstr->IN)[reaidx]);
   ret->put(ret->size(),*(*nstr->OUT)[reaidx]);
   ret->put(ret->size(),*(*nstr->INH)[reaidx]);
   ret->put(ret->size(),*(*nstr->CAT)[reaidx]);
   return ret;
}

/*string* Net::toString(){
	string *st=new string();
	appendToString(*st);
	return st;
}
void Net::appendToString(string &st){
	nstr->appendToString(st);
	reaMask->appendToString(st);
}*/
std::string Net::printToString(){
	size_type i;
   std::string s;
	char charstr[50];
   sprintf(charstr,"Net (%i,%i)\n",(int)nstr->numR,(int)nstr->numC);
	s.append(charstr);
	s.append("IN: [");
	for (i=0;i<nstr->IN->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->IN)[i]->printToString());
	}
	s.append("]\nINST: [");
	for (i=0;i<nstr->INST->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->INST)[i]->printToString());
	}
	s.append("]\nOUT: [");
	for (i=0;i<nstr->OUT->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->OUT)[i]->printToString());
	}
	s.append("]\nOUTST: [");
	for (i=0;i<nstr->OUTST->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->OUTST)[i]->printToString());
	}
	s.append("]\nCAT: [");
   for (i=0;i<nstr->CAT->size();i++){
      sprintf(charstr," %i:",(int)i);
      s.append(charstr);
      s.append((*nstr->CAT)[i]->printToString());
   }
   s.append("]\nINH: [");
	for (i=0;i<nstr->INH->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->INH)[i]->printToString());
	}
	s.append("]\nCMPI: [");
	for (i=0;i<nstr->CMPI->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->CMPI)[i]->printToString());
	}
	s.append("]\nCPRE: [");
	for (i=0;i<nstr->CPRE->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->CPRE)[i]->printToString());
	}
	s.append("]\nCSUC: [");
	for (i=0;i<nstr->CSUC->size();i++){
      sprintf(charstr," %i:",(int)i);
		s.append(charstr);
		s.append((*nstr->CSUC)[i]->printToString());
	}
	s.append("]\nREV: ");
	s.append(nstr->REV->printToString());
	s.append("\nMASK: ");
	s.append(reaMask->printToString());
	return s;
}
/*!
	\brief prints the network to stdout

*/
void Net::print(){
   std::cout << this->printToString() << std::endl;
//	printf("%s\n",this->printToString()->c_str);fflush(stdout);
}
void Net::print_matrix(){
   size_type i,j;
   array<array<std::string> > els;
   std::string el;
   size_type maxlen=0;
   int pos;
   float st;
#ifdef NETUSECOLOR
   bool ansi=true;
#else
   bool ansi=false;
#endif
   if (nstr->numC>60) {
      std::cout << "to large for printing" << std::endl;
      return;
   }
   els.resize(nstr->numC);
   for (i=0;i<nstr->numC;i++){
      for (j=0;j<nstr->numR;j++){
         els[i].resize(nstr->numR);
         el="";
         st=0;
         std::stringstream ss;//create a stringstream
         
         pos=(*nstr->IN)[j]->find(i);
         if (pos>=0) {st-=(*(*nstr->INST)[j])[pos];}
         pos=(*nstr->OUT)[j]->find(i);
         if (pos>=0) {st+=(*(*nstr->OUTST)[j])[pos];}
         ss << st;
         if (ansi) {
            if (st<0) {el+=COL_RED;}
            else if (st>0) {el+=COL_GREEN;}
            else {el+=COL_GREY;}
            st=fabs(st);
         }
         if (st) el+=ss.str();
         if (ansi) el+=COL_RESET;
         if (ansi) el+=COL_BOLD;
         if ((*nstr->INH)[j]->find(i)>=0) el+="i";
         if ((*nstr->CAT)[j]->find(i)>=0) el+="c";
         if ((*nstr->CMPI)[j]->find(i)>=0) el+="I";
         if ((*nstr->PRP)[j]->find(i)>=0) el+="P";
         if (ansi) el+=COL_RESET;
         if (el.length()==0) el="0";
         if (el.length()>maxlen) maxlen=el.length();
         els[i][j]=el;
      }
   }
   for (i=0;i<nstr->numC;i++){
      for (j=0;j<nstr->numR;j++){
         std::cout.width(maxlen);
         std::cout << std::right << els[i][j];
         std::cout << " ";
      }
      std::string nm=(*nstr->names)[i];
      if (nm.length()==0){
         std::stringstream ss;
         ss << "c" << i;
         nm=ss.str();
      }
      std::cout << "| " << nm << std::endl;
   }
}

